Hydrodynamic flow patterns and synchronization of beating cilia 
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We calculate the hydrodynamic flow field generated far from a cilium which is attached to a 
surface and beats periodically. In the case of two beating cilia, hydrodynamic interactions can lead 
to synchronization of the cilia, which are nonlinear oscillators. We present a state diagram where 
synchronized states occur as a function of distance of cilia and the relative orientation of their beat. 
Synchronized states occur with different relative phases. In addition, asynchronous solutions exist. 
Our work could be relevant for the synchronized motion of cilia generating hydrodynamic flows on 
the surface of cells. 
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Many eucaryotic cells possess cilia, which are motile, 
whip- like structures on the cell surface [ij, |2| • While cer- 
tain cells such as sperm swim using a single cilium (in 
this case called flagellum) other cells use several or a large 
number of cilia for propulsion in a fluid. Epithelial cells 
such as those in the respiratory tract use densely cili- 
ated surfaces to transport fluids. Hydrodynamic flows 
generated by cilia also play a key role during the mor- 
phogenesis of higher organisms. In mammals, the left- 
right symmetry of the embryo is broken with the help of 
nodal cilia which generate a hydrodynamic flow to one 
side which transports signaling molecules and breaks the 
symmetry |3j, |J, |5j . All these cilia are based on the same 
structure, called the axoneme, which is built of a very 
regular arrangement of protein filaments called micro- 
tubules in a cylindrical geometry. Motility is achieved 
by the action of a large number of dynein motor proteins 
which generate forces in the cilium while consuming a 
chemical fuel. As a consequence, the cilium produces 
periodic deformations in three dimensional space. 

Cilia in different organisms differ mainly by their 
length and their pattern of beating. They typically mea- 
sure few micrometers up to a few tens of micrometers in 
length and around 150 — 300 nm in diameter. While the 
beat of sperm tails is typically a planar, sometimes heli- 
cal wave, more complex, three dimensional, asymmetric 
beating patterns occur typically in cells which propel flu- 
ids along surfaces. In this case, the working cycle of a 
cilium consists of a fast, upright, effective stroke and a 
recovery stroke which brings the cilium more slowly back 
to the original position on a path closer to the surface, 
see Fig. QJi. Nodal cilia break the left-right symmetry of 
developing embryos by generating rotatory motion in a 
plane tilted with respect to the surface to which they are 
attached 0. 

In this letter, we discuss the hydrodynamic flow field 
generated by a single cilium which beats in an asymmet- 
ric pattern while attached to a surface. While the near 




FIG. 1: a) Beating pattern of a cilium. The effective stroke 
and the recovery stroke are different (arrows) and a net flow 
is generated, b) Simplified representation of the cilium by a 
small sphere, moving on a tilted elliptical trajectory. 



field of the hydrodynamic flow depends on the detailed 
beating pattern and the whip-like geometry of the cil- 
ium, the far field has general features which only depend 
on a few parameters characterizing the symmetry of the 
beat. Using the flow field generated by a beating cilium, 
we study the hydrodynamic interaction between two cilia. 
We discuss conditions that lead to synchronization of the 
two cilia by hydrodynamic interactions. 

Our minimal model of the ciliary beat (Fig. [IJ)) cap- 
tures the essential features and symmetries - the differ- 
ence between effective stroke and recovery stroke of whip- 
like beating as well as the tilted rotatory motion of nodal 
cilia. We replace the cilium by a small sphere of radius a 
(essentially describing the center of mass position of the 
cilium), which moves on a fixed trajectory in the vicinity 
of a planar surface (defined as the plane z = 0). The tra- 
jectory of the bead is elliptic, the phase of the oscillation 
(the position along the trajectory) is described by an an- 
gle 4>. The asymmetry of the ciliary beat is reflected in 
the fact that the two principal axes (denoted A and B) 
are different, and that the ellipse is tilted with respect 
to the surface as described by the parameter C. The 
position of the sphere is given by 
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here, x® and y° describes the position on the surface at 
which the cilium is attached. For more detailed descrip- 
tions of the ciliary beat see, e.g., refs. @jE1- 

The viscous, over damped fluid around the cilium can 
be described by the Stokes equation for the hydrody- 
namic flow field v 



t?V 2 v = VP 



(2) 



where the pressure P is a Lagrange multiplier to impose 
the constraint of incompressibility, V • v = 0. The far- 
field generated by a moving sphere in the vicinity of a 
surface can be calculated by studying the solution of the 
Stokes equation to a delta-distributed force Fj at position 
x i = i x i,Ui: z i) which is of the form v(x) = G(xj,x)Fj. 
The tensor G consists of a stokeslet G s , which describes 
the velocity field around a small isolated particle, and an 
image, required to satisfy the no-slip boundary condition 
on the surface 0. The image is located at the particle 
position, mirrored over the boundary plane, Xj = x, — 
2zi& Zl and consists of an anti-stokeslet, a source-doublet 
G D and a stokes-doublet G SD , so that G(xj, x) = G s (x- 
Xi ) - G s (x - x j: ) + 2z 2 G D (x - Xi) - 2 Zi G SD (x - X;) with 
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The leading behavior of G which describes the far-field 
(for r ^> Zi, z) can be written as 

3 / COs2 ^ S * n ^ C0S ^ ^ 

G(xj,x) » v sin /3 cos 8 sin 2 /3 I (6) 
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where tan 8 = (yi — y)/ — a;) and r 2 = (xi — x) 2 + (j/j — 
y) 2 . For constant height z, G decays as as G cx r~ 3 . 

We assume that the active mechanism in the cilium 
generates a tangential force ft — /jtj. In a simplified 
model, it is a linear function of the velocity Vj = Ujtj 
described by /j = /o — KUj. Here, tj is a normalized vec- 
tor parallel to tj = dxi/d<f>i. Note that the total force 
Fj acting on the sphere is in general not parallel to the 
direction of motion since normal forces arise to estab- 
lish the constraint which keeps the sphere moving on the 
ellipsoidal track. The force Fj thus obeys 



tj • Fj — fo 
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This force is balanced by hydrodynamic friction Fj 
7iVj, where the friction matrix is given by 8] 



7,; = 7(Xj) 




(8) 



Here 70 = 67rrya denotes the Stokes friction of a small 
sphere with radius a and the second term the corrections 
due to proximity of the plane. The force- velocity relation 
0, balanced by the friction force leads to the equation 
of motion for the phase of oscillation 



7(xj)tj + retj • tj 



(9) 



The resulting motion of the sphere is periodic in time. In 
the limit of small radius a, the friction coefficient becomes 
independent of the height over the surface and the period 
is T = 221 ~ (>-+7o)i ^ w here £ denotes the trajectory 
length. The phase <pQ of the sphere as a function of time 
can be written as 



<po{u>t) = ut + K sin(u>t) + L sin(2ojt) + 



(10) 



The coefficient L is related to the eccentricity of the el- 
lipse (if the particle moves at a constant speed, the pa- 
rameter cj) has a variable time derivative) and thus de- 
pends on the geometry of the trajectory. The coefficient 
K oc aC/D 2 , describes variations of the velocity due to 
variations of the friction (JHJ resulting form varying dis- 
tance of the particle from the surface. 

We can now obtain the hydrodynamic far-field of a 

beating cilium v(x, t) = G(xj(t),x)7j(xj(t))vj(t). Av- 
eraging over one period, we define the net flow v(x) = 
J Q T v(x, t)dt = If § G(xj, x)7j(xj)dxj. It only depends 
on the period and the shape of the trajectory, but not 
on details in the phase 4>o(t). Because the sphere causes 
a stronger fluid motion when it is further away from the 
surface (during the working stroke) than when it is closer 
(recovery stroke), it causes a net fluid flow in y direction. 
Examples of flux lines of the time-averaged fluid flow gen- 
erated in this way are displayed in Fig. [3 Far from a 
sphere moving around the center at x® — 0, y° = 0, the 
average velocity field is given by 



9iraBC yz 



(11) 



Note that it only depends on the sphere radius a, the 
period T and the projected area of the trajectory to the 
y-z plane, it BC. Corrections to this flow field can be ne- 
lected in the far field since they decay at least as C(|x| 3 ) 



We now turn to the case of two beating cilia which in- 
teract hydrodynamically. The force F 2 the second sphere 
exerts on the fluid at position x 2 generates hydrodynamic 
flows at position xi and thus influences the force Fi ex- 
erted by the first sphere: 



(12) 



Fi = 71 vi -G(x 2 , Xl )F 2 



This equation can be used to define the dynamics of mo- 
tion of both spheres along their ellipsoidal trajectories 
in the following manner. We assume again that each 
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FIG. 2: a) Flux lines describing the time averaged fluid flow 
generated by a sphere moving on a tilted elliptical path shown 
in the x — y plane for z/D = 2 (other parameters: A/D = 
0.5, B/D = 1, C/D = 0.5, a/D = 0.1). The curves were 
obtained numerically using the flow field based on Eqns. JHJ- 
JSj. b) Same flow field displayed in the y-z plane, for x = 0. 
Black lines correspond to low fluid velocity, yellow lines to 
high velocity. The dashed lines indicate the projections of 
the sphere trajectory. 




FIG. 3: a) Top view of the arrangement of two elliptical tra- 
jectories which represent two beating cilia. The distance r and 
the angle /3 are indicated, b) The state diagram as a function 
of the distance vector x — x\ — x?, V = V2 — Vi, determined 
from the numerical solution of the full model equations. Three 
different regions are indicated. In the region of asynchronous 
beat, two frequencies occur. Two synchronous states can be 
distinguished with equal (Sync + ) and opposite (Sync - ) phase 
of both cilia in the limit of large separation r. 



terized by the phase lag S = (fa ~ <j>x), where the brack- 
ets denote a time average over one period. Regions of 
synchronous states and a region of asynchronous states 
are indicated in the state diagram of Fig. [3Jd. The two 
regions of synchronous states correspond to equal and 
opposite phases in the limit of infinite separation r. The 
values of 8 are displayed in Fig. 0] as a function of the 
angle (3 for different distances r = |xj — x°| between the 
two cilia. For some intervals of (3, no synchronized so- 
lution occurs and the spheres rotate in an asynchronous 
manner with two different frequencies. Note that even 
if both cilia have identical properties, the arrangement 
shown in FigEK breaks the symmetry and both cilia can 
have different preferred frequencies when they interact. 

The existence and stability of synchronous states can 
be studied analytically. First, using Eq. I|12l) . we find for 
large r a 



sphere is driven by a force obeying Eq. (JJJ. Eq. i|12|) 
together with the corresponding equation for F2 then 
uniquely determines the tangential velocities of both 
spheres Vi = tj ■ Vj. This allows us to set up the equa- 
tions of motion that determine the time derivatives fa 
and fa as functions of fa and fa. The problem of two 
coupled phase oscillators is equivalent to the Kuramoto 
model |10j |. a classical model for describing synchroniza- 
tion phenomena |llj |. Performing numerical solutions to 
these dynamic equations, we find that for certain param- 
eter values, after long times the motion of both spheres 
becomes periodic with the same frequency and the oscil- 
lations thus synchronize. A synchronized state is charac- 



/o = ti • (t( x i) + kI)vi - ti • 7(xi)G(x 2 ,xi)7(x 2 )v 2 

(13) 



Here, we have used the fact that for large r S> a, the 
interactions decay as G oc r -3 , while the friction terms 
scale with the sphere size, 7 oc a. Using this relation 
between the tangential velocities, we can study synchro- 
nization. For an isolated sphere, the phase as a function 
of time is denoted (f>o(t) lfTU|) . We denote by <j)\{t) the 
phase of a sphere 1 when interacting with sphere 2. The 
variation AT of the oscillation period T of sphere 1 can 
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FIG. 4: The phase difference S between two spheres on elliptic 
trajectories in the stable stationary state as a function of their 
relative position in polar coordinates (angle j3, distance r) ob- 
tained by numerical solutions to the full dynamic equations. 
In the limit r — > oo analytic far-field calculations show that 
only equal and opposite phases (5 = and ±7T, respectively) 
occur. For smaller r, nontrivial angles 5 occur. Synchronous 
states become unstable where the lines end (black dots). Be- 
tween these points motion is asynchronous. 



then be written as 
AT = 
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ti • 7(xi)G(x 2 , xi)7(x 2 )t 2 v 2 dt (14) 



where the last expression becomes correct in the limit 
of large r/a. In this limit, the variation AT can be 
calculated using the Ansatz Xi = x° + x(<po(u>t)) and 
x 2 = X2+x(0o(^^+^)) resulting in a function AT(r, j3, 8). 
A synchronized steady state exists if the change in pe- 
riod due to interactions is the same for both spheres, 
i.e., AT(r,(3,8) — AT(r, f3 + ir, — 8), where we have taken 
into account that exchanging both spheres corresponds 
to /3 — ► (3 + 7T and 8 — > —8, see Fig. Further- 
more, a synchronized steady state is locally stable if 
d/d8[AT(r,/3,S) - AT(r,/3 + tt, -8)] < 0. In order to 
find explicit expressions for AT, we perform a system- 
atic expansion in the small parameter K. In the limit of 
vanishing radius a, the friction matrix 7 becomes inde- 
pendent of the height z. The velocity of a single sphere 
then becomes constant along the trajectory (K = 0). As 
a consequence, the corresponding variation ATo of the 
period, resulting from Eq. (|14l> . becomes symmetric with 
respect to 8 — > — 8 as well as with respect to (3 — > f3 + tt 
and the condition that both particles have the same pe- 
riod is satisfied trivially, regardless of the phase lag 8. In 
order to find nontrivial synchronized states, we have to 
go to first order in K where 



2(7o + K)£ 2 r 3 

(AT S (8) = AT s (—8) represents symmetric terms that are 
irrelevant for synchronization). At this order in a and r, 
the condition for a stable synchronized steady state is 
fulfilled with 8 = for < < tt/2 and tt < < 3tt/2 
and with 8 = n for 7r/2 < /3 < n and 3ir/2 < f3 < 2n. 
This is consistent with our numerical solutions in the 
limit of large r, see Fig. 0] For smaller distance r, the 
contribution of terms proportional r -4 becomes relevant. 

To conclude, we have shown that the hydrodynamic 
far held around a periodically beating cilium has generic 
properties which do not depend on the detailed beat- 
ing pattern. Its main features can thus be generated by 
a sphere moving on a tilted elliptical trajectory. This 
captures the asymmetry of the ciliary beat. As a con- 
sequence, the resulting flow is a time periodic pattern 
superimposed on a time-independent net flow. If the 
beating cilium is perturbed by external forces, this in- 
terferes with the internal force generating process and 
affects the instantaneous angular velocity of oscillations. 
We capture this effect by a linear relationship between 
external force and the velocity v of the sphere. The force 
exerted on one cilium by the hydrodynamic flow gener- 
ated by the other couples both oscillators. The result- 
ing synchronization phenomena depend on the distance 
r between cilia but also the parameters a/D and a/C 
which characterize the difference of effective and recov- 
ery strokes. No synchronization occurs if the motion is 
rotationally symmetric or helical [12T Il3| . 

A natural extension of our study is the generalization 
to a periodic lattice of cilia attached to a surface. Such 
situations correspond to some epithelial cells which gen- 
erate fluid flows on their surface and to microorganisms 
such as Paramecium. In all these cases, the ciliary beat is 
organized in metachronal waves which result from hydro- 
dynamic and steric interactions between cilia. In different 
cells, these waves of ciliary strokes propagate along the 
surface in different directions relative to the one defined 
by the working stroke [l4l[T5| . Other examples for hydro- 
dynamically stirred surfaces are realized in experiments 
where bacteria are attached to a surface by their flag- 
ella which are driven by rotatory motors [lg ■ Each bac- 
terium in this system could be represented by a rotating 
sphere as in our description. In principle, the synchro- 
nization effects between two rotating elements discussed 
here could be studied in such experiments. The extension 
of our study to many hydrodynamically coupled cilia and 
the description of metachronal waves and other collective 
modes is left for future work. 
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